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ABSTRACT 

We present steady solutions of high speed viscous flows over blunt bodies using a multidomain 
Chebyshev spectral collocation method. The region within the shock layer is divided into subdo- 
mains so that internal layers can be well- resolved. In the interiors of the subdomains, the solution 
is approximated by Chebyshev collocation. At interfaces between subdomains, the advective terms 
are upwinded and the viscous terms are treated by a penalty method. The method is applied to 
five flows the Mach number range 5-25 and Reynolds number range 2,000 - 83,000, based on nose 
radius. Results are compared to experimental data and to a finite difference result. 
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Introduction 


Spectral methods are now applied to a wide variety of fluid dynamics problems^. They are 
characterized by global, high order orthogonal polynomial approximations. For smooth enough solutions they 
give exponential order convergence. The high order approximation of the advective terms has low dissipation 
errors, which is important in viscous flow calculations. 

Compared to finite difference methods, fewer grid points are needed by spectral methods to obtain 
highly accurate solutions. As an example, Pruett and Streett^ compared the results of a Chebyshev method for 
the solution of the compressible boundary layer equations to that of a second order finite difference method. To 
compute the wall temperature for an adiabatic wall, they found that 77 points were required for the spectral 
method to obtain the same five digit accuracy for which the finite difference method required 400. 

On the negative side, the use of global polynomial approximations makes it difficult to solve problems 
with multiple length scales. These scales might be associated with the size of an aerodynamic body, the 
distance between shock waves, the size of entropy and vorticity layers, and the size of viscous boundary layers. 
A standard spectral collocation method would require a complex mapping and a large number of grid points to 
solve such problems. Refinement of the grid to resolve local features would have to be handled by the 
mappings. The resulting problem would be numerically very stiff. Very small time-steps and large numbers of 
iterations would be required to converge to steady-state. 

To ameliorate the difficulties of standard spectral methods, multidomain spectral methods have been 
developed to solve incompressible viscous (see Ref. 1 for a discussion) and, now, compressible-^ flows. 
Multidomain methods still give the exponential error convergence characteristic of the spectral methods, but 
also allow for local refinement of the grid. 

In this paper, we apply to the solution of the steady blunt body problem a new multidomain Chebyshev 
collocation method developed for the solution of the compressible Navier-Stokes equations^. The important 
feature of the method is the treatment of the interfaces between subdomains. At a subdomain interface, the 
advective terms in the equations are treated by an upwind characteristic correction procedure developed for the 
in viscid Euler gas-dynamics equations. ^ The viscous terms are approximated by a weighted average of the 
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terms on either side plus a penalty term proportional to the jump in the viscous fluxes. The result is a method 
with communication limited to the passing of boundary data only. 

The method is applied to five examples. The first four examples are of sphere, cylinder and blunt 
cone geometries at Mach number 5.73 and 10.6. For these solutions, the surface pressure or the heat transfer are 
compared to experiments. The final problem is the solution of a Mach 25 flow over a hyperbolic cone. This 
problem was solved earlier by a single domain spectral method.^ The results are compared to a finite difference 
Thin Layer Navier-Stokes calculation and to the single domain solution. 

The Blunt Body Problem 

In this paper we solve for viscous supersonic flow about a blunt nose cone at zero angle of attack. Fig. 
1 shows a diagram of the situation. Spectral approximations of this problem using single domains have been 
performed by Kopriva^ and Yasuhara, Nakamura and WangW Because of the symmetry along the centerline 
of the body, we need to solve only the upper half of the full flow problem. For high Mach and Reynolds 
numbers, the bow shock can be viewed as a discontinuity 7 By fitting the shock as a boundary to the flow, then, 
we consider the solution only in the shock-layer. 
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Fig. 1 Diagram of the blunt body geometry. The gridded area represents the computation 
region subdivided into multiple subdomains 


We approximate the compressible Navier-Stokes equations in the non-conservative form 
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The parameter % = 1 for axisymmetric flow and £ = 0 otherwise. 

In these equations, we have scaled the dimensional pressure, density, velocity and temperature ip*, p*, 
q *, T*) by 


P=P*/Pr'f 


P = P*!Pref 

q = (m»v) = q * /-{p^f / p r< / (3a) 
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where 9t is the gas constant We scale the dimensional space and time by 


x = x*IL (3b) 

* = 'sjprtf t P rtf * * IL> 

where L represents the length scale. For the blunt body problems computed here, we choose p re f and p T ef to be 
the free stream values and I to be the nose radius. 

The temperature is computed by the perfect gas equation of state, which bee (Hues 


T=plp (4) 

under the scaling (3). The dimensionless variables in (1) are the standard Reynolds Number (Re), which is 
scaled to the characteristic length, the Prandtl Number (Pr) and the Mach Number (Af). 

To compute the viscosity coefficient, we use the Sutherland law 


H* 1 + 198.6/ 
n rff r+198.6/7^ 


(5) 


where T re f is measured in °R. We assume that the Prandtl number is constant, so the same formula is used to 
compute the dimensionless heat conductivity, k. 


The Multidomain Spectral Method 

The multidomain method that we use is presented in Kopriva^. As before, we consider non- 
overlapping and “conforming” approximations only. By conforming we mean that subdomains intersect only 
along a side or at a point, and that grid lines must be continuous between subdomains. However, in this paper, 
the geometry is allowed to change with time (because of the moving shock) and axisymmetric geometry is 
considered. The moving grid requires the interface conditions to change at each time step. 

The domain decomposition for the blunt body problem is performed as described in Kopriva^ We 
subdivide the region between the shock and the body into multiple non-overlapping subdomains, which 

collectively cover the region between the shock and the body. (See Fig. 1) Each subdomain is bounded by four 
curves. In the direction normal to the body, r, the bounds are r^(s 9 t) and Time dependent 

subdomain boundary positions are included to allow the shock boundary to move with time. For convenience 


only, we require the subdomains to be bounded in the body tangential direction, s y by straight lines normal to the 
body surface. 

Given the boundary curves on a subdomain, a mapping from ( r,s,t ) ( X t Y,T) is used: 

r-rj inCM) 
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For each subdomain the (6) transformation causes the equations (1) to be rewritten as 
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In (7), the advective terms are 
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For convenience, we have left off the superscripts that denote the subdomain. The viscous terms are 
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Here, the Jacobian, J, is 


J - x x y r x r y x 

The viscous fluxes in the mapped coordinate are defined to be 

K = y r F t - XyG, 

G t =-y x F t + x x G' 

for q = u, v, orp , and 

F u = n(4 (X x u x + Y x u r )~ 2(X y v x + f A )) /3-^ 
G u = F y = li{X x v x + Y x v Y + X y u x + Y y u r j 
G y = n(4(X,v x + Y x v r )-2(X y u x + Y y u r ))l3-^- 
F P = *{XJ x + Y x T y ) 

G p = K(X y T x + Y y T r ) 
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Finally, 
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We do not multiply through by the Jacobian of the transformation to create new dependent variables, since J 
can be discontinuous across interfaces between subdomains. 

On each subdomain, G*, we place a grid of points 




i = 0,1,2,..., N k 
)) j = 0 ,\, 2 ,...,M k 


(13) 


and assign gridpoint values of the solution unknowns, e.g. p* •, etc. Derivatives of a solution variable, q (= u,v, 
or p), are then approximated from the grid point values, qu, by 
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( 14 ) 


(A,- i« 

J=0 
/=0 

where d i ■ is the Chebyshev derivative matrix as written in Ref. 1, pg. 69. We use matrix multiplication rather 


than FFT’s to compute the derivative approximation. FFT’s are more efficient for large numbers of points. 
However, the number of points used here is usually chosen to be fewer than 32, which is around the break-even 
number of points where the FFT’s become faster. 

The derivative approximations (14) are substituted into eqns. (8) - (12) to give a system of ordinary 
differential equations, which we denote as 


and, for q = each of p , u, or v, 



dql 

dt 


= -(< + (< 


(15a) 


(15b) 


These equations are integrated in time with a second order Runge-Kutta rule. 

At interfaces between the subdomains, conditions are applied to allow waves to propagate through 
them, while weakly enforcing continuity of the viscous flux. To describe the interface method, we consider two 
subdomains as shown in Fig. 2, one on the left (L) and the other on the right (R) of the interface. 
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Fig. 2. Diagram of a two subdomain decomposition 


At each interface point, we replace the spatial derivatives in (7) by their discrete Chebyshev 
approximations, (14). Two equations for each grid point on the interface are generated, depending on whether 
derivatives are computed from the left or from the right. For the density, the equation is purely advective, and 
we write 




(16a) 


For q = each of p, u, v, the equations on either side of the interface are written so that the viscosity is 
approximated by a weighted average of the viscous terms computed on either side of the interface. The 
continuity required of the viscous flux is enforced weakly by a penalty proportional to the jump in the flux: 


% - -KC, * “MX ♦ 


^ - -Ml + «MX *f,Ml * Me - ej 


(16b) 


The weights, a, /J, S depend on the relative sizes of the subdomains, and on the numbers of points in the 
subdomains: 
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The quantities, w* = 




, are the Clenshaw-Curtis quadrature weights deflned 


on the interval [0,1]. Thus, the coefficient of the penalty term grows as 0(N^), and in the limi t as N (16b) 
enforces continuity of the flux. To update the solutions, the equations (16) are integrated with a second order 
Runge-Kutta rule. 
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An approximation similar to (16) can be written for comer points found at the intersection of four 
subdomains^. At such points, the weighted average is taken over the four values of the diffusion terms, and 

penalties are applied for the jumps in the diffusive fluxes across each face. 

The results of the integration in time are preliminary values of the solutions, p, p, w, V on either side 


of the interface. The preliminary values have used one sided, but not upwind, approximations for the advective 
terms at the interface points. Thus, these solutions must be corrected to account for wave propagation across the 
interface. From w,v, preliminary values of the normal and tangential velocities, vv 1 and vv 2 , are formed. The 

preliminary values are then corrected to have the proper advective domain of dependence in the manner 
described in Ref. 3. For instance, if the flow is subsonic and from left to right across the interface, the choice 
of four bicharacteristics of the hyperbolic portion of the equations leads to the corrected interface values 



The superscript n refers to the linearization to the previous time level. The velocities u and v can then be 
computed from the corrected normal and tangential velocities, Wj and 

The boundary conditions and their approximations are described in detail in Ref. 9. We will only 
review the basic ideas here. To compute the wall pressure, we use a compatibility equation that relates the 
pressure and the normal velocity along the in viscid characteristic direction. Thus, the pressure is required to 
satisfy a combination of both the energy and momentum equations at the body surface. If the surface is 
adiabatic, this equation uses the normal temperature gradient set equal to zero when the second derivative is 
computed. The density is computed from the pressure and the temperature either from the equation of state, for 
an isothermal surface, or, for an adiabatic surface, by a spectral extrapolation that enforces the zero normal 
temperature gradient. 
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The shock boundary is treated as a discontinuity by shock fitting since at high Mach and Reynolds 
numbers, the shock can be considered a sharp discontinuity. The shock-fitting procedure is described in detail in 
Ref. 4 for the inviscid case. The sharp shock assumption implies that the viscous terms are negligible near the 
shock, so we ignore them in the calculation of the downstream pressure. 

At the outflow boundary, we impose the viscous conditions by setting the normal viscous fluxes to be 
zero when the second derivatives are computed. Since the streamwise flux is small in high Reynolds number 
flow, this has little effect on the solution. In the subsonic portion of the boundary layer, we also specify the 
pressure to be a fixed value, and typically use the inviscid pressure. That pressure can be obtained from a 
Newtonian flow assumption, which is reasonably accurate for hypersonic flows. It can also be obtained from 
the solution of the inviscid equations, which we use to start the viscous computation. 

Unlike the solutions presented in Ref. 10, no artificial viscosity or spectral filtering are used with this 
approximation. Filtering of the solution either directly or by the addition of artificial viscous terms is 
undesirable for three reasons. First, the overall accuracy of the approximation is decreased, since the number of 
accurate frequencies is reduced. Second, high frequency information is altered. Finally artificial viscous effects 
can be introduced in the solution. In the absence of shocks or other discontinuities in the computed region, 
which we avoid here by shock fitting, the need for explicit filtering of smooth solutions often indicates the use 
of poor boundary procedures that introduce high frequency errors into the solution. These errors destroy both 
the solution accuracy and, in steady-state problems, the convergence rate. 
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Solutions 


Five problems are solved using the method described above. A summary of the flow parameters is 
shown below in Table I. The problems include both two-dimensional and axisymmetric flows at a range of 
Reynolds and Mach Numbers. The solutions for problems I-IV are compared to experimental data. The last 
problem compares the computed results with a finite difference solution^. 


Table I. Parameters for the Flow Test Cases. 


Case 

Geometry 

Moo 

Re 

Too 

T\y/Too 

Purpose 

I 

Sphere 

10.6 


85.21R 

Adiab. 

Compute C D 

II 

Cylinder 

5.73 

2050 

71.37R 

71.37 

Heat Transfer 

III 

15° Sphere/Cone 

10.6 

8.33X10 4 

85.21R 

Adiab. 

C D 

IV 

15° Sphere/Cone 

10.6 

1.65xl0 5 

85.21 R 

6.22 

Heat Transfer 

V 

5° Hyperboloid 

25 


486 R 

25 

Compare with Thin Layer Solutior 


Case I. 

The first problem that we solve is the M = 10.6 flow over a sphere at Reynolds number 83.000 based 
on nose radius, using a grid with eight subdomains. Fig. 3 shows the computed contours of the pressure, Mach 
number and temperature along with the grid used. The narrow strip of subdomains along the surface of the 
sphere was chosen to better resolve the boundary layer. A coarse grid is all that is needed to resolve the 
essentially inviscid flow outside of the boundary layer. Convergence behavior of the maximum residual over all 
flow variables, plotted in Fig. 4, shows that steady-state is obtained in double precision (64 bit) arithmetic. 

The solution contours are smooth and pass continuously through the interfaces. Normal direction 
profiles for the pressure, velocity and temperature and their derivatives are shown in Fig. 5, at s = 0.78, 
approximately halfway along the surface of the sphere. 
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Mach 



Fig. 3. Solution contours and grid for Case I. 
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Finally, Fig. 6 plots the pressure coefficient 


P^-P. 

z-pja 

2 y 


along the surface of the sphere measured as horizontal distance from the leading edge (See Fig, 1). The result is 
compared to the inviscid solution computed in Ref. 5, and to the experimental data of Cleary^. The viscous 
computation clearly approximates better the data than the inviscid computation. 
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Normal Distance (r) 

Fig. 5 Solution and derivative profiles at s = 0. 78 for Case I. 
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d Iql/dr dT/dr 










Fig. 6 Pressure coefficient for Case l as a function of horizontal distance from the 

leading edge. 

Case II. 

For the second problem, the heat transfer over a cylinder placed in a M = 5.73 free stream was 
computed, again using eight subdomains. This problem was also solved using a single domain in Ref 9. Fig. 7 
shows the contours of the pressure, Mach number and the temperature along with the grid. As before, the 
solutions are smooth through the interfaces. The heat transfer along the surface of the cylinder is shown in Fig. 
8. The results are compared to the experiments of Tiewfik and Geidt^ and to the finite difference results of 
Gnoffo^. Smoothness of the solution is implied by the continuity of the normal derivative of the temperature 
across the interfaces. 
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Fig. 7 Solution contours and grid for Case II. 
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Case III. 

A more complex situation is shown in Fig. 9, which depicts the Mach 10.6 flow over a spherically 
blunted 15° half-angle cone. The solution contours in Fig. 9 show both the growth of the boundary layer along 
the cone surface and the shrinking of the entropy layer generated by the shock. The solution was computed on a 
grid of 14 subdomains. Fig. 10 shows a comparison of the pressure coefficient along the surface of the cone to 
the inviscid calculation^ and to the experimental results^. Marked by arrows on Fig. 10 are the positions of 
the subdomain interfaces. 
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Pressure 



Fig. 9 Solution contours and grid for Case 111. 
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Fig. 10 Pressure coefficient along body surface for Case III 

Case IV. 

Case IV is similar to Case III except that the temperature is kept fixed along the surface of the cone. 
Fig. 1 1 shows the heat flux normalized to the value of the nose. The computed solution is compared to the 
experimental results^. 




Fig. 11. Heat flux for Case IV. 
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CaseV 


As a final example, we make more detailed comparisons of the results of a M = 25 flow over a 5° 
hyperbolic cone to those of the shock-fitted finite difference code of Ref 12. The finite difference calculation 
used 50 points in each direction. Since the ideal gas law, (4), is used, the calculation is meant only for a 
comparison of the numerical results. This calculation used eight subdomains. Subdomains along the body used 
13 points in the normal direction, while those along the shock used seven. 

The calculation presented here required 140,000 time steps (starting with the inviscid solution) for the 
density to reduce the maximum residual in all flow variables 12 orders of magnitude. (See Fig. 12) At 0.285 sec 
per time step, this means that the calculation required approximately 11 hrs of CPU time on an IBM RS/6000 
320H. For comparison purposes, this superscalar machine runs at approximately 4-5 megaflops. 



0 50 100 150 

n/1000 

Fig. 12. Convergence of the maximum residual for Case V. 


20 



We will first compare the shock position and the wall values of the pressure coefficient, normalized 
heat transfer and skin friction as a function of the distance along the cone. We will also compare the profiles of 
the flow quantities in the direction normal to the cone surface at four stations along the length of the cone. We 
begin by comparing the shock position, as shown in Fig. 13. Plotted is the normal distance of the shock from 
the body as a function of arc length along the body. Also compared to the multidomain solution is a single 
domain solution ^ computed on a 20x20 grid so that the number of points in the normal direction is the same. 

0.8 

0 ) 0.7 

o 

5 0.6 

g 0.5 

% 0.4 

T3 

| 0-3 

35 0.2 

0.1 

0 1 2 3 4 5 6 

s 

Fig. 13 Shock standoff distance comparing multidomain spectral, 
single domain spectral and finite difference solutions for Case V . 

Fig. 14 compares the pressure coefficient, the heat flux Q> and the skin friction along the body surface 
as a function of the distance from the nose. 

Variations of the pressure, temperature and velocity in the direction normal to the body surface are 
compared at four stations along the body in Fig. 15. The four stations correspond to the symmetry line and 
three interfaces that intersect the body. In all cases, we see excellent agreement between the spectral solution 
and the finite difference solution, and smooth variations through the interfaces. 
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Conclusions 


We have used a multidomain spectral collocation method to solve the compressible Navier-Stokes 
equations for two-dimensional and axisymmetric supersonic flows over blunt bodies. A main advantage of the 
method over a single domain method is the ability to distribute easily both subdomain boundary locations and 
the number of points in each subdomain. Thus, features such as boundary layers can be easily resolved. 

At interfaces between subdomains, the advection terms are treated with the characteristic correction 
method used in Ref. 3 for the inviscid Euler equations. In this way, waves propagate through the interfaces in 
much the same way as they would for an inviscid flow. The diffusion terms are treated with a penalty method 
The penalty weakly enforces the continuity of the viscous flux. 

The method was applied to five problems. The computations show that the method is stable and that it 
converged in all cases to a smooth steady-state solution without the need for artificial viscosity or filtering. The 
results were compared to experimental data, to a finite difference solution, and to a single domain spectral 
method. 

Spectral methods are a high resolution methods and for smooth enough solutions give exponential 
order accuracy. Phase and dissipation errors are exponentially small. Along with the high accuracy, however, is 
the fact that the pointwise work is significantly higher than that required by low order finite difference methods. 
This means that the computations are relatively expensive. Hanley^ has presented evidence for one space 
dimension that if high enough accuracy is required, a spectral method can be as much as a factor of seven more 
efficient than a second order MacCormack scheme. Nevertheless, there is still a serious need for research on 
methods to accelerate the convergence rate of spectral methods. At this time, we consider the method described 
here to be appropriate for nose solutions. For longer bodies, the nose solution could be used as a starting 
solution for a marching code, since the approximation is defined by its Chebyshev interpolant even between 
grid points. 
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